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ABSTRACT 


The  need  to  compute  eigenvalues  and  eigenvectors 
for  large  order  symmetric  matrices  arises  frequently 
in  the  finite  element  approach  to  static  structural 
analysis  as  well  as  in  the  applied  work  of  many  other 
fields  in  the  Navy's  varied  research  program.  GIVENS 
is  an  out -of -core  FORTRAN  callable  eigensystem  subroutine 
which  computes  all  the  eigenvalues  and  specified  eigen- 
vectors for  a large  order  real  symmetric  matrix  using 
modifications  of  the  GIVENS  tridiagonal izat ion  procedure, 
the  QR  method  of  Francis,  and  inverse  iteration.  GIVENS 
has  been  put  together  from  subroutines  extracted  from 
the  NASTRAN  (NAsa  STRuctural  Analysis)  program.  The  use 
of  GIVENS  presupposes  that  the  user  can  attach  a special 
FORTRAN  compiler  and  a certain  subroutine  library 
presently  available  on  the  CDC  6000  series  computers  at 
DTNSRDC. 


INTRODUCTION 

One  of  the  long  range  projects  of  the  Computation,  Mathematics,  and 
Logistics  Department  has  been  the  development  of  mathematical  subroutines 
suitable  for  use  in  the  computer-aided  structural  analysis  of  ships. 

Many  unrelated  efforts  in  both  government  and  industry  have  resulted  in 
computer  programs  that  treat  particular  classes  of  structural  problems. 
These  programs  often  involve  the  solution  of  similar  mathematical  problems 
but,  since  the  solutions  are  reached  independently,  the  efficiency  and 
accuracy  of  the  various  algorithms  used  may  vary  greatly.  The  need  to 
coordinate  these  diverse  efforts,  to  develop  improved  methods  of  more 
general  applicability,  and  to  produce  more  comprehensive  programs  for 
solving  Navy  structural  problems  became  obvious.  A project  was  therefore 
established  to  coordinate  research  efforts  involving  mathematical  and 


1 


computational  methods  in  the  area  of  structural  mechanics  and  to  integrate 
the  work  of  mathematicians , computer  specialists,  and  structural  engineers 
in  this  field. 

The  present  considerable  interest  in  the  finite  element  approach  to 
structural  analysis  is  evidenced  by  the  widespread  use  of  NASTRAN  (NAsa 
STRuctural  ANalysis  program) ^ and  other  such  programs.  Fundamental  to  the 
finite  element  approach  to  vibration  analysis  is  the  solution  of  the 
matrix  equation 

(K  - AM)U  = 0 

for  natural  frequencies  (A)  and  normal  modes  (U) . If  both  K and  M are 
real  and  symmetric  and  M is  positive  definite,  this  problem  can  be  reduced 
to  one  of  computing  the  eigensystem  of  a real  symmetric  matrix  C.  However, 
the  order  of  the  problem  is  often  so  large  that,  even  when  advantage  is 
taken  of  sparsity,  it  is  not  feasible,  and  often  not  even  possible,  to 
store  the  matrix  C in  the  core  (high  speed  direct  access)  memory  of  a 
computer,  let  alone  to  compute  C's  eigensystem  in  core.  The  need  to 
compute  eigensystems  for  such  large  order  real  symmetric  matrices  also 
occurs  quite  frequently  in  the  applied  work  of  other  scientific  disciplines 
in  the  Navy's  far-ranging  research  program. 

The  NASTRAN  program  has  mathematical  "modules"  to  cope  with  this  and 
other  mathematical  or  computational  problems  which  arise  in  the  course  of 
the  NASTRAN  structural  analysis.  To  be  sure,  these  modules  can  be  accessed 
directly  for  a mathematical  problem  as  such  but  then  DMI  (Direct  Matrix 
Input)  cards  are  required  for  the  matrices.  Unfortunately,  this  DMI  card 
procedure  is  not  realistic  for  large  order  matrices.  Accordingly,  an 
investigation  was  undertaken  to  see  whether  some  of  the  mathematical 
capabilities  of  these  modules  could  be  extracted  from  the  NASTRAN  program 
and  recast  in  the  form  of  FORTRAN  callable  subroutines  with  a more  efficient 
input  facility  for  large  order  matrices.  The  present  GIVENS  subroutine  was 
put  together  in  this  way  with  the  assistance  of  Mr.  Myles  Hurwitz  (Code 
1844)  who  removed  the  appropriate  subroutines  from  the  NASTRAN  READ  (Real 

' "The  NASTRAN  User's  Manual  (bevel  15.0),"  NASA  SP-222,  Scientific  and 
Technical  Information  Office,  National  Aeronautics  and  Space  Administration, 
Washington,  D.C.,  .June  1972. 


Eigenvalue  Analysis  for  the  Displacement  method)  module  along  with  the 
required  general  library  input  and  output  subroutines  trom  NASTRAN. 


THEORY 

This  section  briefly  describes  the  procedure  used  by  GIVENS  to  com- 
pute the  eigensystem  of  a real  symmetric  matrix.  For  a more  detailed 
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presentation  the  reader  is  referred  to  the  NASTRAN  Theoretical  Manual  or 
to  the  authors  cited  below.  The  GIVENS  procedure  is  as  follows: 

(1)  The  matrix  is  reduced  to  a real  symmetric  tridiagonal  matrix  using 

3 4 

Wilkinson's  modification  of  the  GIVENS  method. 

(2)  All  the  eigenvalues  of  this  tridiagonal  matrix  are  computed  using 
Ortega  and  Kaiser's  modification"’  of  the  QR  method  of  Francis.^  Since  the 
two  matrices  are  similar,  these  are  also  the  eigenvalues  of  the  original 
matrix. 

(3)  Certain  eigenvectors  [those  corresponding  to  the  smallest  eigenvalues 
of  those  eigenvalues  in  a certain  range)  of  the  tridiagonal  matrix  are 

7 

computed  using  a method  developed  by  Wilkinson.  The  Gram-Schmidt  algo- 

g 

rithm  is  used  to  insure  eigenvector  orthonormality  in  the  case  of 
repeated  eigenvalues. 

2 "The  NASTRAN  Theoretical  Manual  (Level  15.0),"  NASA  SP-221(01), 
Scientific  and  Technical  Information  Office,  National  Aeronautics  and  Space 
Administration,  Washington,  D.C. , December  1972,  pp.  10.1-1,  10.2-15. 

Wilkinson,  J.H.,  "The  Algebraic  Eigenvalue  Problem,"  Oxford  University 
Press,  1965,  pp.  506-510. 

4 

Givens,  W. , "Numerical  Computation  of  the  Characteristic  Values  of  a 
Real  Symmetric  Matrix,"  Oak  Ridge  National  Lab.,  ORNL-1574,  1954. 

5 T 

Ortega,  J.M.  and  H.F.  Kaiser,  "The  LL1  and  QR  Methods  for  Symmetric  Tri- 

diagonal  Matrices,"  Computer  J. , Vol.  6,  No.  1 (Jan.  1963),  pp.  99-101. 

^ Francis,  J.G.F.,  "The  QR  Transformation,  a Unitary  Analogue  to  the  LR 
Transformation,"  Computer  J.,  Vol.  4,  No.  3 (Oct.  1961)  and  No.  4 
(Jan.  1962). 

Wilkinson,  J.H. , "The  Calculation  of  the  Eigenvectors  of  Codiagonal 
Matrices,"  The  Computer  J. , Vol.  1,  1958,  pp.  90. 

g 

Bodcwig,  F,.,  ’"Matrix  Calculus,”  Interscience  Publisher,  1959. 
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(4)  The  eigenvectors  of  the  original  matrix  are  obtained  from  those  of 

the  tridiagonal  matrix  by  multiplying  each  tridiagonal  eigenvector  by  the 
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series  of  orthogonal  rotation  matrices  determined  by  GIVENS'  method. 


THE  GIVENS  SUBROUTINE 

The  calling  sequence  for  GIVENS  is 

CALL  GIVENS (N ,NV , JTIME , LFREC ,HFREC) 

The  value  of  the  first  integer  argument  N is  the  order  n of  the  matrix. 

Now  n must  be  entered  as  the  value  of  some  variable  N rather  than  as  a 
numerical  constant,  since  the  argument  N is  redefined  in  GIVENS.  This 
restriction  does  not  apply  to  the  other  arguments  of  GIVENS. 

The  values  of  the  second  integer  parameter  NV  and  of  the  fourth  and 
fifth  real  arguments  LFREC  and  HFREC  are  assigned  as  follows:  If  the  user 

wishes  to  compute  the  eigenvectors  of  the  k smallest  eigenvalues,  the 
value  of  NV  is  k and  the  values  of  LFREC  and  HFREC  are  immaterial.  If  the 
user  wishes  to  compute  the  eigenvectors  corresponding  to  eigenvalues  X in 
the  interval  a < X < b,  the  value  of  NV  is  zero  and  the  values  of  LFREC 
and  HFREC  are  a and  b,  respectively.  The  number  of  eigenvectors  computed 
is  returned  through  the  argument  NV.  If  the  user  merely  wishes  all  the 
eigenvalues  and  no  eigenvectors,  the  value  of  NV  is  zero  and  the  common 
value  of  LFREC  and  HFREC  is  c,  where  c is  not  an  eigenvalue.  (The  reader 
is  reminded  that  GIVENS  always  computes  all  the  eigenvalues.) 

The  value  of  the  third  integer  argument  .JTIME  is  a time  estimate 
(integral  number  of  seconds)  for  the  eigensystem  computation.  One  such 
estimate  is  given  by  the  formula 

JTIME  = N3  + lSxlO'5  . 

GIVENS  reads  the  rows  (or  columns)  of  the  matrix  from  tape  8,  a row 
(or  a column)  at  a time.  GIVENS  writes  the  eigensystem  output  on  tape  9 
in  the  following  way.  The  first  record  contains  the  N eigenvalues  and 
the  succeeding  NV  records  contain  the  specified  eigenvectors  (if  any). 


At  present  GIVENS  is  available  in  the  form  of  a pre -compiled 
permanent  file  on  the  CDC  6400.  The  user's  calling  program  must  be 
compiled  with  the  FORTRAN  compiler  originally  used  to  compile  NASTRAN  to 
ensure  compatibility  with  certain  NASTRAN  subroutines.  GIVENS  also  requires 
certain  subroutines  from  the  SYSMISC  library.  The  sample  setup  below 
illustrates  the  control  cards  required  for  GIVENS  on  the  CDC  6400. 

JOB  CARD 

CHARGE  CARD 

ATTACH, SYSMISC 

L I BRARY , SYSMISC, SYSIO. 

ATTACH , F I LE , G I VENSE I GNSYM , I D=CAD6 . 

ATTACH ,FTN , FTN3POLV340 , ID=CSYS . 

FTN  CARD 

RFL  (to  field  length  on  job  card) 

LOAD, FILE. 

LGO. 

END  OF  RECORD  CARD 

PROGRAM  MAIN  ( 1 NT! IT  .OUTPUT  ,TAPES= INPUT  ,TAPE6=OUTPtJT , 

TAPE 8 , TAPE9) 

REAL  LFREC 


N = 

CALL  G EVENS (N,NV ,JTIME, LFREC, HFRF.C) 

The  field  length  CM  for  the  job  card  may  be  determined  as  follows: 

CM^,  the  amount  of  core  memory  required  by  GIVENS , is  given  by  the  formula 

CM6  = N2  + 6N  + 2655 

where  CM^  is  expressed  as  an  octal  number.  Let  CM^  be  an  estimate  for  the 
field  length  of  the  calling  program  for  GIVENS.  Then  the  total  field 
length  CM  is  the  octal  sum 

CM  = CNV,  - CM6 

(CM^  can  be  determined  more  exactly  from  the  load  map  of  a successful  run.) 

Similarly  the  time  T (integral  number  of  seconds)  for  the  job  card  is 
obtained  from  the  decimal  sum 

T = Tm  + JTIME 

where  T,,  is  a time  estimate  for  the  calling  program  and  .JTIME  is  the 
third  argument  of  GIVENS. 
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Let  A be  the  real  symmetric  matrix  of  order  n = 10,  120,  10  with  10 

on  the  diagonal  and  1 elsewhere.  A is  a positive  definite  matrix  with 

n 9 

n-1  repeated  eigenvalues  of  9 and  one  eigenvalue  of  9+n.  Let  Bn  be  the 
real  symmetric  banded  matrix  of  order  n = 10,  120,  10  obtained  by 
imposing  a bandwidth  of  10  on  is  likewise  positive  definite  but 

with  distinct  eigenvalues. 

Table  1 contains  the  times  (in  seconds)  for  tridiagonal izat ion  and 
the  tridiagonal  eigenvalue  computation  for  A^.  Table  2 contains  the  same 
data  for  B . Four  subroutine  implementations  of  the  basic  algorithm 
outlined  previously  are  compared  in  these  tables:  GIVENS;  TRI4  and 

IMQL2;10  EHOUSS  and  EQRT25;11  TRED1  and  IWTQL1 . 12  The  times  were  obtained 
on  the  CDC  6400  at  DTNSRDC.  The  GIVENS  times  were  obtained  with  the 
Scope  3.3  operating  system,  and  the  other  times  were  obtained  with  the 
Scope  3.4  operating  system. 


Westlake,  J. , "A  Handbook  of  Numerical  Matrix  Inversion  and  Solution 
of  Linear  Equations,"  John  Wiley  § Sons,  Inc.,  1968,  p.  141. 

10  Nikolai,  P.J.  and  N.-K.  Tsao,  "The  ARL  LINEAR  ALGEBRA  LIBRARY  Hand- 
book," Interim  Report,  1 September  1973-31  March  1974,  ARL  TR  74-0106 
(July  1974),  Aerospace  Research  Laboratories,  Air  Force  Systems  Command, 
Wright -Patterson  AFB. 

11  The  IMSL  LIBRARY  3 Reference  Manual,  Vol.  I,  Edition  4 (FORTRAN  2.4) 
CDC  6200/6400/6500/6600/7600  (1975),  revised  November  1974,  International 
Mathematical  and  Statistical  Libraries,  Inc.,  Houston,  Texas. 

12 

Smith,  B.T.,  et  al , ''Matrix  Eigensystem  Routines  - EISPACK  Guide," 
Lecture  notes  in  Computer  Science,  Vol.  6,  Springer-Verlag,  New  York 
(1974) . 
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OBSERVATIONS 

The  GIVENS  program  was  undertaken  with  a twofold  purpose:  first,  tc 

see  how  difficult  it  is  to  extract  a NASTRAN  module,  recast  it  as  a 
stand-alone  FORTRAN  callable  subroutine,  and  then  evaluate  this  subroutine 
as  such;  and  second,  to  provide  an  out -of -core  eigensystem  capability  for 
users  of  our  computer.  This  work  has  indicated  that  it  is  not  too 
difficult  to  remove  a NASTRAN  module  and  make  a subroutine  out  of  it  -- 
assuming  familiarity  with  the  NASTRAN  program.  In  the  present  instance 
the  NASTRAN  expertise  was  provided  by  Mr.  Myles  Hurwitz  of  Code  1844. 
Module  removal  should  not  be  attempted  without  a NASTRAN  consultant . 

No  out-of-core  symmetric  matrix  eigensystem  solvers  wrere  at  hand  to 
compare  with  GIVHNS.  Accordingly,  suitable  in-core  subroutines  were 
chosen  from  the  various  subroutine  libraries  for  comparison  purposes.  The 
times  for  these  comparison  runs  are  presented  in  Tables  1 and  2.  /Ml 
four  subroutines  produced  virtually  identical  results  as  regards  accuracy. 
It  should  also  be  borne  in  mind  that  GIVENS  requires  significant  amounts 
of  "PP"  time.  The  behavior  of  the  ARL  subroutine  TRI4  is  surprising  and 
warrants  investigation. 

GIVENS  is  an  experimental  program.  Because  it  was  extracted  from 
NASTRAN,  it  is  perhaps  a bit  more  complicated  than  if  it  had  been  written 
directly.  As  expected,  GIVENS  has  a time  disadvantage  when  in-core 
computation  is  possible.  GIVENS  is  meant  for  sparse  matrices  too  large 
to  fit  in  core. 
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SUflROUT  INE  GIVt NS ( N,NV ,JT IME, L FRtC , HFR.C) 

INTEGER  T,T1,T2,T3 

INTEGER  SYSBUF 

REAL  LFREC 

EXTERNAL  REAO, WRITE 

DIMENSION  MCBC7) ,IZ(1>  ,1ITuE<1) 

DIMENSION  N AM ( 2 ) 

EQUIVALENCE  <Z,IZ> , (TITLE, II TlC) 

COMMON  Z(l> 

COMMON  /PACKX/  IN, I OUT ,11 , NN, I NCR 
COMMON  /GIVN/  TITLE  ( lb  C ) 

COMMON  /ST  I ME/  I T I ME 

COMMON  /UNPAKX / 1 1 OUT , II 1 , NNN,  1 1 NCR 
COMMON  /SYSTEM/  SYSBUF, NOUT 
COMMON  /MS&X/  NMSG 
DATA  NAM/uHGIV£,**MNS  / 


SET  J3  END  OF  OPEN  CORE  FOR  GINO  INOICE3 

NFILE  S = ( N J.  CF  GINO  FlLtS*l)  *(N3RMST+’i)=21*b3=1.32J 

NFILES=1 323 
LC  ORE  = KORSZ ( Z) 

LCORt^LSORE+NFILES 
J=  LDORE-NFILES  *1 
DO  20  K=J,lCORw 
2 8 I Z ( X ) = 0 

SET  UP  AOQRESS  763  FOR  THOSE  ROUTINES  USING  GORjZ,  NOT  KORSZ 
LCORF=P 

CALL  FI£LDLN(LCORE) 

LC0RE  = LC0R'--NFILES 
CA._  XSTORr (LCORt, 1,763) 

CALL  XSTORl  (lCORE,  1 ,633) 


IT  1ME  = JT  I Me. 

IIT_E(13H=N 
TITLE ( 1 ] 2 ) = L FR - C 
TITLL(i:E)=hFREC 
IITLt (1 J7) =NV 
TITLE (11L)=C 
MCR(  1)  =3li*. 

MS  3( 2) -z 
MCB(  3)  =N 
MCB(-) =6 
MS  B( 5 ) =1 
MCB(6) =3 
MCB( 7) =0 
NN  - N 

LLCORE=KORSZ(Z) 

IS  TO  RE  =MA  XO  ( 19  *N,  7 *N  «-2*SY$9U-  ) 
IMOPE=ISTORl-LLCORE 

12 


UK n,*T ■ ~~r. T t ** 
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IF  < IMORE)  27,2  7. 3o 
?’  NZ  = L LC  ORE 
AN  = N 

flZ  = SZ-F“N-2*SVSoU>" 

M = SOBT  ( az> 

A<sAN-AM*1. 

AMB-15. E-6 

fi»i<=mc  .E-6 

Tl=9.0*AN*AN,AN*AMn*ll.0*AN*AN«APAK 
T2  = 0.  3 

IF  (AK.LE.0.0)  GO  TO  1101 

T2  = APAKMAK*(AK*1.0)*(AN*2.J>-<  AN*  AM*  3.  J)*AK*(AK*1.3>*.i> 

1 f A<*  ( AKH.  Q)  M AK*A<*1.  J> /b.  Z ) 

1131  CONTINUE 
AV  = NV 

T3  = b«C*AN*AN*AV*AiB*3.3*AV*AV*AN*AM3*6.*AV*AN*A:)AK*AV*AV*£DAK 

T=T1*T2*T3 

n=  an 

M=  AM 

WRITE (NOUT, 1000) T ,N,M 
CALL  TMTOGO(I) 

IF (I .GE. T ) GO  TO  101 

I P 1 = - 3 0 

IFILE=T 

CALL  MES  AGE  (IPi,IFILE»NAM) 

1C1  CONTINUE 

10(10  FORMAT (60H0***  USER  I N rO  - ^ A T I 0 N M^bSAS.  2(Jla,  GIVENS  TIME  ESTlMATt. 

• IS  ,IU,  9H  SECONDS.  / 

• 3&X ,lbMPROBLEM  SIZE  IS  ,14, 

* 5UH,  SPILL  WILL  OCCUR  FOR  THIS  CORE  AT  A PROBLEM  SIZE  OF  , 

* 13, 2H  .) 

l:ore=korsz(z> 

IBUF=LCORE-SYSBUF*l 
LCO»E=IBUF-l 
IMORE=N-LCORE 
IF  ( N-LCORE) 29,23,30 
2B  CALL  G0PENC3C4,Z<IBUF)  ,1) 

REWIND  3 

DO  5 J=1,N 

READ!  8)  (Z(  I)  ,I  = 1,N> 

CALL  PACK(Z(1) ,304, WRITE, MCB) 
b CONTINUE 
REWIND  8 

CALL  CLOSE( 304,1) 

CALL  W RT  T RL (MC  3 ) 

CALL  VALVEC 
NV  = 1 1 TLE  (107) 

Ir  ( NV.LE. 1 ) GO  TO  62 

IFiN.EQ. l)GO  TO  62 

m;b(  i ) = 2 ; 3 

MCP ( 2) =N 

MO  R ( 3 ) - N 

M0B(  4)  =3 

MC  8 ( 5 ) = 2 

M0B(  6)  =2 
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MCB(7) =1U0UG./N 
CALL  WRTTRL(MC3) 

MAA=2J3 
EPSII  = 1.  t-*» 

MC8(  1) =3U6 
CALL  KDTRKMCd) 

CALL  REAO4(201,MCb,3i.6,EPSII,MAA> 

62  CONTINUE 
I T 0 J T = 1 
1 1 1 = 1 
NNN=N 

ii n:r=i 

CALL  GOPtN(2Ul,Z(  IGUF)  ,U> 

CALL  REA0<2C1,Z<1)  ,N,0  ,M)  , RETURNS (40 ,45) 

CALL  CL0SE(2ul,l> 

RE  NINO  9 

WRITE<9>  (7 < L K)  ,LK=1,N> 

IF  (NV)  34 , 3 5,  34 
3-  CALL  GOPEN(3u5,Z(IBUF>  ,3) 

DO  1C  J=1,NV 

CALL  UNPACK (31):?  ,Z ( 1)  ,READ)  ,RETURNS<15> 

GO  TO  ?5 
ie  DO  2C  1 = 1, N 

2 . z ( i > = : . 

2b  WRITE (9)  <Z(LK>  , L K = 1 , N ) 

1.  CONTINUE 
REWIND  9 

CALL  CLOSE < 3C5 ,1) 
lb  IF  (NHSG.GT.O)  CALL  MSGHRT 
RE  TURN 

3:  WRITE ( 6,31)IM0R£,IM0RE 

31  FORMAT  (1H1/1*»(1HC/)  /37X,19MlN SUFFICIENT  CJRE.,i3,3M  / ,0a, 
1 >♦  1 h MORE  DECIMAL  / OCTAL  LOCATIONS  rlQU1R_D.) 

ST  0° 

- . WRITE!  6,ui) 

-1  FORMAT ( 1H1/1^( lHj/) /61X,i2HE0F  ON  lAMA.) 

STOP 

4 E WR  I T E ( 6 , *♦  b ) 

»F  FORMAT  (lHl/l4(lhG/)/61X,l2HEJR  ON  LAMA.) 

STOP 

END 
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